R script of Data-driven Regularized Adaptive Risk Forecasting

pkg_list = c('zoo', 'tseries', 'fGarch','PEIP','tidyverse','gridExtra', 'gdata', 'xtable','vioplot')
# ensure existing required packages are up to date:
#update.packages(ask=FALSE, oldPkgs=pkg_list)
# Install packages if needed
for (pkg in pkg_list)
{
  # Try loading the library.
  if ( ! library(pkg, logical.return=TRUE, character.only=TRUE) )
    {
         # If the library cannot be loaded, install it; then load.
        install.packages(pkg, repos = "http://cran.us.r-project.org")
        library(pkg, character.only=TRUE)
  }
}

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Loading required package: timeDate

## Loading required package: timeSeries

## 
## Attaching package: 'timeSeries'

## The following object is masked from 'package:zoo':
## 
##     time<-

## Loading required package: fBasics

## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   1.0.0
## ✓ tidyr   1.1.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks timeSeries::filter(), stats::filter()
## x dplyr::lag()    masks timeSeries::lag(), stats::lag()

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.

## 

## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.

## 
## Attaching package: 'gdata'

## The following object is masked from 'package:gridExtra':
## 
##     combine

## The following objects are masked from 'package:dplyr':
## 
##     combine, first, last

## The following object is masked from 'package:purrr':
## 
##     keep

## The following object is masked from 'package:stats':
## 
##     nobs

## The following object is masked from 'package:utils':
## 
##     object.size

## The following object is masked from 'package:base':
## 
##     startsWith

## 
## Attaching package: 'xtable'

## The following object is masked from 'package:timeSeries':
## 
##     align

## The following object is masked from 'package:timeDate':
## 
##     align

## Installing package into '/Users/hoqueme/Library/R/4.0/library'
## (as 'lib' is unspecified)

## also installing the dependency 'sm'

## 
## The downloaded binary packages are in
##  /var/folders/j9/dfgsrh6n401466rqzvbmlmq80000gn/T//RtmpYIN9pl/downloaded_packages

## Loading required package: sm

## Package 'sm', version 2.2-5.6: type help(sm) for summary information

#library(fBasics)
rm(list=ls())# Remove objects from enviornment
dateStart <- "2005-01-01"               
dateEnd <- "2018-11-28"

# Apple data
APPLE <- get.hist.quote(instrument="AAPL",start = dateStart, end=dateEnd, quote = c("AdjClose"),
                        retclass="zoo")

## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.

## time series starts 2005-01-03
## time series ends   2018-11-27

# Microsoft data
MSFT <- get.hist.quote(instrument="msft",start = dateStart,end=dateEnd,quote = c("AdjClose"),
                       retclass="zoo")

## time series starts 2005-01-03
## time series ends   2018-11-27

# IBM data
IBM <- get.hist.quote(instrument="ibm",start = dateStart,end=dateEnd,quote = c("AdjClose"),
                      retclass="zoo")

## time series starts 2005-01-03
## time series ends   2018-11-27

# S&P 500 index, which is one of the leading stock market indices for US equity 
SP500 <- get.hist.quote(instrument="^gspc",start = dateStart, end=dateEnd,quote = c("AdjClose"),
                        retclass="zoo")

## time series starts 2005-01-03
## time series ends   2018-11-27

#Google data
GOOG <- get.hist.quote(instrument="GOOG",start = dateStart, end=dateEnd,quote = c("AdjClose"),
                        retclass="zoo")

## time series starts 2005-01-03
## time series ends   2018-11-27

data <- merge(SP500,IBM,MSFT,APPLE,GOOG) # closing price data
return.cc = diff(log(data))
ret.cc<-as.data.frame(tail(return.cc, 3000)) # latest 3000 observations
p <- seq(0.001, 0.01, 0.0001) #p for VaR
main.names<-c("SP500", "IBM", "MSFT","AAPL","GOOG")
#sample correlation
rho.cal<-function(X){
  rho.hat<-cor(sign(X-mean(X)), X-mean(X))
  return(rho.hat)
}
ret.0.cc<-ret.cc
for(j in 1:5){
  ret.0.cc[, j]<-ret.cc[, j]-mean(ret.cc[, j])
}
rho<-apply(as.matrix(ret.0.cc), MARGIN=2, FUN=rho.cal)
#calculate degree of freedom
nu<-rep(0, 5) 
for(i in 1:5){
  fun <- function (x) rho[i]*(x-1)*beta(x/2,1/2)-2*sqrt(x-2)
  nu[i] <- uniroot(fun, c(2, 8))$root
}
number<-ncol(data)
acf.s<-rep(0, number)
acf.abs<-rep(0, number)
acf.sq<-rep(0, number)
for(j in 1:number){
  acf.s[j]<-acf(ret.0.cc[, j], plot=FALSE)$acf[2]
  acf.abs[j]<-acf(abs(ret.0.cc[, j]), plot=FALSE)$acf[2]
  acf.sq[j]<-acf(ret.0.cc[, j]^2, plot=FALSE)$acf[2]
}
corr<-data.frame(apply(ret.cc, 2, mean), apply(ret.cc, 2, sd), 
                 apply(ret.cc, 2, kurtosis),apply(ret.cc, 2, skewness),acf.s, acf.abs, 
                 acf.sq, rho, nu)
rownames(corr)<-main.names
colnames(corr)<-c("mean", "sd","kurtosis","skewness","series", "abs", "sq", "sign-rho", "df")
xtable(corr, digits=4)

## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:18 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrrr}
##   \hline
##  & mean & sd & kurtosis & skewness & series & abs & sq & sign-rho & df \\ 
##   \hline
## SP500 & 0.0002 & 0.0124 & 10.8749 & -0.3711 & -0.1002 & 0.2861 & 0.2152 & 0.6351 & 2.9877 \\ 
##   IBM & 0.0002 & 0.0139 & 5.9892 & -0.2393 & -0.0075 & 0.2125 & 0.1483 & 0.6894 & 3.6302 \\ 
##   MSFT & 0.0005 & 0.0172 & 9.4608 & 0.1642 & -0.0687 & 0.2374 & 0.1571 & 0.6721 & 3.3673 \\ 
##   AAPL & 0.0009 & 0.0199 & 7.2219 & -0.4408 & 0.0077 & 0.2487 & 0.1832 & 0.6920 & 3.6775 \\ 
##   GOOG & 0.0005 & 0.0180 & 11.6231 & 0.5192 & -0.0115 & 0.1730 & 0.0734 & 0.6622 & 3.2461 \\ 
##    \hline
## \end{tabular}
## \end{table}

ret.0<-tail(ret.0.cc, 1000) #latest 1000 observations
##### normal garch
VaR_rateG=matrix(0, nrow=length(p), ncol=ncol(ret.0))
ES_rateG=matrix(0, nrow=length(p), ncol=ncol(ret.0))
sigma_rate=c()
omega=c()
alpha=c()
beta=c()
for (j in 1:ncol(ret.0))
{
  g = garchFit(~garch(1,1),ret.0[,j],cond.dist="norm",include.mean=FALSE,
               trace=FALSE)
  omega[j] = g@fit$matcoef[1,1]
  alpha[j] = g@fit$matcoef[2,1]
  beta[j] = g@fit$matcoef[3,1]
  sigma_rate[j] =  omega[j] + alpha[j] * ret.0.cc[,j][T]^2 + beta[j] * g@h.t[T]
  for(i in 1:length(p)){
  VaR_rateG[i,j] = -sqrt(sigma_rate[j]) * qnorm(p[i]) * 1000
  ES_rateG[i,j] = -sqrt(sigma_rate[j])*integrate(function(q){q*dnorm(q)},-Inf,qnorm(p[i]))$value/(p[i]) * 1000
  }
}
####### t garch
VaR_ratetG1=matrix(0, nrow=length(p), ncol=ncol(ret.0))
sigma_rate=c()
omega=c()
alphaG=c()
beta=c()
ES_ratetG1=matrix(0, nrow=length(p), ncol=ncol(ret.0))
dfG=c()
for (j in 1:ncol(ret.0))
{
  g = garchFit(~garch(1,1),ret.0[,j],cond.dist="std",include.mean=FALSE,
               trace=FALSE)
  omega[j] = g@fit$matcoef[1,1]
  alphaG[j] = g@fit$matcoef[2,1]
  beta[j] = g@fit$matcoef[3,1]
  dfG[j]=g@fit$matcoef[4,1]
  sigma_rate[j] =  omega[j] + alpha[j] * ret.0[,j][T]^2 + beta[j] * g@h.t[T]
  for(i in 1:length(p)){
    VaR_ratetG1[i,j] =(-1)* qstd (p[i], mean = 0, sd = sqrt(sigma_rate[j]), nu = dfG[j])*1000  
  # using fGarch package to get df
  ES_ratetG1[i,j] = sqrt(sigma_rate[j])*sqrt((dfG[j]-2)/dfG[j])*dt(qt (p[i],dfG[j]), dfG[j])/p[i]*(dfG[j] + (qt (p[i],dfG[j]))^2)/(dfG[j]-1)*1000}
}

ret.0.cc<-as.data.frame(ret.cc)
for(j in 1:5){
  ret.0.cc[, j]<-ret.cc[, j]-mean(ret.cc[,j])
}

alpha.choose<-function(ret.cen, alpha=seq(0.01, 0.3, 0.01), cut.t){
  rho<-rho.cal(ret.cen)
  vol<-abs(ret.cen)/rho
  t<-length(ret.cen)
  MSE_alpha<-rep(0, length(alpha))
  for(a in 1:length(alpha)){
    s<-mean(vol[1:cut.t])
    error<-rep(0, t)
    for(i in 1:t){
      error[i]<-vol[i]-s
      s<-alpha[a]*vol[i]+(1-alpha[a])*s
    }
    MSE_alpha[a]<-mean(error[-(1:cut.t)]^2)
  }
  rmse<-sqrt(min(MSE_alpha))
  alpha.opt<-alpha[which.min(MSE_alpha)]
  return(c(rmse, alpha.opt))
}

alpha.opt<-rep(0, number)
rmse.usual<-rep(0, number)
for(num in 1:number){
  result<-alpha.choose(ret.0.cc[, num], seq(0.01, 0.3, 0.01), 2000)
  alpha.opt[num]<-result[2]
  rmse.usual[num]<-result[1]
}
alpha.opt;rmse.usual

## [1] 0.14 0.05 0.07 0.07 0.06

## [1] 0.008803848 0.013319630 0.015985200 0.014808302 0.015813455

rmse.elastic.cal<-function(ret.cen, alpha, omega, lambda, cut.t){
  rho<-rho.cal(ret.cen)
  vol<-abs(ret.cen)/rho
  #vol.sub<-vol-sd(ret.cen)
  s<-mean(vol[1:cut.t])
  sd.r<-sd(ret.cen)
  if(abs(s-sd.r)>=omega*lambda)  
    s.elastic<- sign(s-sd.r)*(abs(s-sd.r)-omega*lambda)/(1+(1-omega)*lambda)+sd.r
   else s.elastic<-sd.r
  t<-length(ret.cen)
  error<-rep(0, t)
  for(i in 1:t){
    error[i]<-vol[i]-s.elastic
    s<-alpha*vol[i]+(1-alpha)*s
    if(abs(s-sd.r)>=omega*lambda)  
      s.elastic<-sign(s-sd.r)*(abs(s-sd.r)-omega*lambda)/(1+(1-omega)*lambda)+sd.r
    else s.elastic<-sd.r
  }
  rmse<-sqrt(mean(error[-(1:cut.t)]^2))
  return(rmse)
}

lambda<-seq(0, 0.004, 0.0002)
omega<-seq(0, 1, 0.1)
rmse.elastic.opt<-rep(0, number)
omega.opt<-rep(0, number)
lambda.opt<-rep(0, number)
for(num in 1:number){
  rmse.elastic<-matrix(0, nrow=length(omega), ncol=length(lambda))
  for(a in 1:length(omega)){
    for(l in 1:length(lambda)){
      rmse.elastic[a, l]<-rmse.elastic.cal(ret.0.cc[, num], alpha.opt[num], 
                                           omega[a], lambda[l], 2000)
    }
  }
  index.o<-which.min(apply(rmse.elastic, 1, min))
  index.l<-which.min(apply(rmse.elastic, 2, min))
  omega.opt[num]<-omega[index.o]
  lambda.opt[num]<-lambda[index.l]
  rmse.elastic.opt[num]<-rmse.elastic[index.o, index.l]
}
omega.opt; lambda.opt

## [1] 0.2 0.2 0.6 0.3 0.3

## [1] 0.0034 0.0040 0.0022 0.0034 0.0036

rmse.elastic.opt

## [1] 0.008775514 0.013299324 0.015928461 0.014782626 0.015784693

vol.forecast<-function(ret.cen, alpha, omega, lambda, cut.t){
  rho<-rho.cal(ret.cen)
  vol<-abs(ret.cen)/rho
  sd.r<-sd(ret.cen)
  t<-length(ret.cen)
  res<-rep(0, t)
  res.elastic<-rep(0, t)
  s<-mean(vol[(1:cut.t)])
  s.elastic<-sign(s-sd.r)*max((abs(s-sd.r)-
               omega*lambda), 0)/(1+(1-omega)*lambda)+sd.r
  fore<-rep(0, t)
  fore.elastic<-rep(0, t)
  for(i in 1:t){
    res[i]<-ret.cen[i]/s
    res.elastic[i]<-ret.cen[i]/s.elastic
    fore[i]<-s
    fore.elastic[i]<-s.elastic
    s<-alpha*vol[i]+(1-alpha)*s
    s.elastic<-sign(s-sd.r)*
              max((abs(s-sd.r)-omega*lambda), 0)/(1+(1-omega)*lambda)+sd.r
  }
  fore.vol<-s
  fore.vol.elastic<-s.elastic
  return(list(fore.vol, fore.vol.elastic, res, res.elastic, fore, fore.elastic))
}

forecast.vol<-rep(0, number)
forecast.vol.elastic<-rep(0, number)
res<-ret.0.cc
res.elastic<-ret.0.cc
fore.vary<-ret.0.cc
fore.elastic.vary<-ret.0.cc
for(num in 1:number){
  result<-vol.forecast(ret.0.cc[, num], alpha.opt[num], omega.opt[num], 
                       lambda.opt[num], 2000)
  forecast.vol[num]<-result[[1]]
  forecast.vol.elastic[num]<-result[[2]]
  res[, num]<-result[[3]]
  res.elastic[, num]<-result[[4]]
  fore.vary[, num]<-result[[5]]
  fore.elastic.vary[, num]<-result[[6]]
}
forecast.vol; forecast.vol.elastic

## [1] 0.01462323 0.01935135 0.02464491 0.03018339 0.02447984

## [1] 0.01393917 0.01853662 0.02331948 0.02914136 0.02338628

vol.obs<-ret.0.cc
for(num in 1:number){
  vol.obs[, num]<-abs(ret.0.cc[, num])/rho[num]
}
for(num in 1:number){
  plot(1:1000, tail(vol.obs[, num], 1000), xlab="t", ylab="vol", main=main.names[num], lwd=3,
       col="yellow", type="l")
  lines(1:1000, tail(fore.vary[, num], 1000), col="blue", lwd=5)
  lines(1:1000, tail(fore.elastic.vary[, num], 1000), col="purple")
  abline(h=sd(ret.0.cc[, num]), lty=3, col="red", lwd=3)
  legend("top",legend=c("vol", "NP", "Elastic", "sd"), col=c("yellow", "blue", "purple", "red"),
         lty=c(1, 1, 1, 3), lwd=c(3, 5, 1, 3), cex=0.7)
}

for(num in 1:number){
  plot(1:1000, tail(fore.vary[, num], 1000), col="green", lwd=5, type="l")
  lines(1:1000, tail(fore.elastic.vary[, num], 1000), col="purple", lwd=4)
  abline(h=sd(ret.0.cc[, num]), lty=3, col="red", lwd=3)
  legend("top",legend=c( "NP", "Elastic", "sd"), col=c( "green", "purple", "red"),
         lty=c(1, 1, 3), lwd=c( 5, 1, 3), cex=0.7)
}

vol.CI<-matrix(0, nrow=number, ncol=2)
vol.CI.elastic<-matrix(0, nrow=number, ncol=2)
for(i in 1:number){
  vol.CI[i, ]<-c(forecast.vol[i]+qnorm(0.05, 0, 1)*rmse.usual[i]/sqrt(1000),
                 forecast.vol[i]+qnorm(0.95, 0, 1)*rmse.usual[i]/sqrt(1000))
  vol.CI.elastic[i, ]<-c(forecast.vol.elastic[i]+qnorm(0.05, 0,1)*rmse.elastic[i]/sqrt(1000),
                       forecast.vol.elastic[i]+qnorm(0.95, 0, 1)*rmse.elastic[i]/sqrt(1000))
}

scale.p<-function(nu, p) sqrt((nu-2)/nu)*dt(qt (p,nu), nu)/p*(nu + (qt (p,nu))^2)/(nu-1)

VaR.dd<-matrix(0, nrow=length(p), ncol=number)
ES.dd<-matrix(0, nrow=length(p), ncol=number)
VaR.dd.elastic<-matrix(0, nrow=length(p), ncol=number)
ES.dd.elastic<-matrix(0, nrow=length(p), ncol=number)
for(num in 1:number){
  for(i in 1:length(p)){
    VaR.dd[i, num]<-forecast.vol[num]*qstd(p=p[i], nu=nu[num])*(-1000)
    VaR.dd.elastic[i, num]<-forecast.vol.elastic[num]*qstd(p=0.01, nu=nu[num])*(-1000)
    ES.dd[i, num]<-forecast.vol[num]*scale.p(nu[num], p[i])*1000
    ES.dd.elastic[i, num]<-forecast.vol.elastic[num]*scale.p(nu[num], p[i])*1000
  }
}
VaR.dd[91,]; VaR.dd.elastic[91, ]

## [1] 38.30383 51.44188 65.44769 80.22095 64.86404

## [1] 36.51201 49.27607 61.92785 77.45145 61.96643

ES.dd[91, ]

## [1]  59.18422  73.77755  96.32483 114.57789  96.79783

no.penalty<-data.frame(alpha.opt, rmse.usual, forecast.vol,(t(VaR.dd)[,91]), (t(ES.dd)[,91]) )
colnames(no.penalty)<-c("alpha.opt", "RMSE", "Volatility", "VaR", "ES")
rownames(no.penalty)<-main.names
xtable(no.penalty, digits=6)

## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
##   \hline
##  & alpha.opt & RMSE & Volatility & VaR & ES \\ 
##   \hline
## SP500 & 0.140000 & 0.008804 & 0.014623 & 38.303826 & 59.184225 \\ 
##   IBM & 0.050000 & 0.013320 & 0.019351 & 51.441877 & 73.777549 \\ 
##   MSFT & 0.070000 & 0.015985 & 0.024645 & 65.447693 & 96.324834 \\ 
##   AAPL & 0.070000 & 0.014808 & 0.030183 & 80.220948 & 114.577890 \\ 
##   GOOG & 0.060000 & 0.015813 & 0.024480 & 64.864037 & 96.797826 \\ 
##    \hline
## \end{tabular}
## \end{table}

penalty<-data.frame(omega.opt, lambda.opt,rmse.elastic.opt, forecast.vol.elastic,
                    t(VaR.dd.elastic)[, 91],(t(ES.dd.elastic)[,91]))
colnames(penalty)<-c("omega.opt", "lambda.opt","RMSE", "Volatility", "VaR", "ES")
rownames(penalty)<-main.names
xtable(penalty, digits=6)

## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrr}
##   \hline
##  & omega.opt & lambda.opt & RMSE & Volatility & VaR & ES \\ 
##   \hline
## SP500 & 0.200000 & 0.003400 & 0.008776 & 0.013939 & 36.512009 & 56.415641 \\ 
##   IBM & 0.200000 & 0.004000 & 0.013299 & 0.018537 & 49.276072 & 70.671367 \\ 
##   MSFT & 0.600000 & 0.002200 & 0.015928 & 0.023319 & 61.927854 & 91.144393 \\ 
##   AAPL & 0.300000 & 0.003400 & 0.014783 & 0.029141 & 77.451450 & 110.622275 \\ 
##   GOOG & 0.300000 & 0.003600 & 0.015785 & 0.023386 & 61.966433 & 92.473676 \\ 
##    \hline
## \end{tabular}
## \end{table}

VaR.dd.CI.m<-matrix(0, nrow=number, ncol=4)
ES.dd.CI.m<-matrix(0, nrow=number, ncol=4)
for(num in 1:number){
  VaR.dd.CI.m[num, 1:2]<-(-1000)*qstd(p=0.01, nu=nu[num])*vol.CI[num,]
  VaR.dd.CI.m[num, 3:4]<-(-1000)*qstd(p=0.01, nu=nu[num])*vol.CI.elastic[num, ]
  ES.dd.CI.m[num, 1:2]<-(1000)*scale.p(nu[num], 0.01)*vol.CI[num, ]
  ES.dd.CI.m[num, 3:4]<-(1000)*scale.p(nu[num], 0.01)*vol.CI.elastic[num, ]
}

rho.res<-apply(res, 2, rho.cal)
rho.res.elastic<-apply(res.elastic, 2, rho.cal)
nu.res<-rho.res
nu.res.elastic<-rho.res.elastic
for(num in 1:number){
  fun <- function (x) rho.res[num]*(x-1)*beta(x/2,1/2)-2*sqrt(x-2)
  nu.res[num] <- uniroot(fun, c(2, 8))$root
  fun <- function (x) rho.res.elastic[num]*(x-1)*beta(x/2,1/2)-2*sqrt(x-2)
  nu.res.elastic[num] <- uniroot(fun, c(2, 8))$root
}
nu.res.elastic

## Adjusted.SP500   Adjusted.IBM  Adjusted.MSFT Adjusted.APPLE  Adjusted.GOOG 
##       4.649606       3.927487       4.045986       4.644023       3.728411

nu.res

## Adjusted.SP500   Adjusted.IBM  Adjusted.MSFT Adjusted.APPLE  Adjusted.GOOG 
##       4.323821       3.765093       3.859418       4.472556       3.652099

VaR.dd.res<-matrix(0, nrow=length(p), ncol=number)
ES.dd.res<-matrix(0, nrow=length(p), ncol=number)
VaR.dd.elastic.res<-matrix(0, nrow=length(p), ncol=number)
ES.dd.elastic.res<-matrix(0, nrow=length(p), ncol=number)
for(num in 1:number){
  for(i in 1:length(p)){
    VaR.dd.res[i, num]<-forecast.vol[num]*qstd(p=p[i], nu=nu.res[num])*(-1000)
    VaR.dd.elastic.res[i, num]<-forecast.vol.elastic[num]*qstd(p=p[i], nu=nu.res.elastic[num])*(-1000)
    ES.dd.res[i, num]<-forecast.vol[num]*scale.p(nu.res[num], p[i])*1000
    ES.dd.elastic.res[i, num]<-forecast.vol.elastic[num]*scale.p(nu.res.elastic[num], p[i])*1000
  }
}
nu.res

## Adjusted.SP500   Adjusted.IBM  Adjusted.MSFT Adjusted.APPLE  Adjusted.GOOG 
##       4.323821       3.765093       3.859418       4.472556       3.652099

VaR.dd.elastic.res[91,]

## [1] 36.55100 49.15726 61.74651 76.42124 62.13669

ES.dd.elastic.res[91,]

## [1]  49.08741  68.84120  85.76310 102.65919  88.36805

VaR.dd.CI.res<-matrix(0, nrow=number, ncol=4)
ES.dd.CI.res<-matrix(0, nrow=number, ncol=4)
for(num in 1:number){
  VaR.dd.CI.res[num, 1:2]<-(-1000)*qstd(p=0.01, nu=nu.res[num])*vol.CI[num,]
  VaR.dd.CI.res[num, 3:4]<-(-1000)*qstd(p=0.01, nu=nu.res.elastic[num])*vol.CI.elastic[num, ]
  ES.dd.CI.res[num, 1:2]<-(1000)*scale.p(nu.res[num], 0.01)*vol.CI[num, ]
  ES.dd.CI.res[num, 3:4]<-(1000)*scale.p(nu.res.elastic[num], 0.01)*vol.CI.elastic[num, ]
}
VaR.CI<-data.frame(VaR.dd.CI.m[, 1:2], VaR.dd.CI.res[, 1:2],
                   VaR.dd.CI.res[, 3:4], VaR.dd.CI.res[, 3:4])
xtable(VaR.CI)

## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
##   \hline
##  & X1 & X2 & X1.1 & X2.1 & X1.2 & X2.2 & X1.3 & X2.3 \\ 
##   \hline
## 1 & 37.10 & 39.50 & 37.35 & 39.76 & 34.39 & 38.71 & 34.39 & 38.71 \\ 
##   2 & 49.60 & 53.28 & 49.56 & 53.24 & 46.98 & 51.34 & 46.98 & 51.34 \\ 
##   3 & 63.24 & 67.66 & 63.20 & 67.61 & 59.57 & 63.92 & 59.57 & 63.92 \\ 
##   4 & 78.17 & 82.27 & 77.36 & 81.41 & 74.26 & 78.58 & 74.26 & 78.58 \\ 
##   5 & 62.68 & 67.04 & 62.88 & 67.26 & 59.95 & 64.32 & 59.95 & 64.32 \\ 
##    \hline
## \end{tabular}
## \end{table}

ES.CI<-data.frame(ES.dd.CI.m[, 1:2], ES.dd.CI.res[, 1:2],
                   ES.dd.CI.res[, 3:4], ES.dd.CI.res[, 3:4])
rownames(ES.CI)<-main.names
xtable(ES.CI)

## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
##   \hline
##  & X1 & X2 & X1.1 & X2.1 & X1.2 & X2.2 & X1.3 & X2.3 \\ 
##   \hline
## SP500 & 57.33 & 61.04 & 51.00 & 54.30 & 46.19 & 51.98 & 46.19 & 51.98 \\ 
##   IBM & 71.14 & 76.42 & 70.27 & 75.49 & 65.79 & 71.90 & 65.79 & 71.90 \\ 
##   MSFT & 93.08 & 99.57 & 88.95 & 95.16 & 82.74 & 88.79 & 82.74 & 88.79 \\ 
##   AAPL & 111.65 & 117.50 & 104.80 & 110.29 & 99.76 & 105.56 & 99.76 & 105.56 \\ 
##   GOOG & 93.55 & 100.05 & 90.01 & 96.27 & 85.26 & 91.48 & 85.26 & 91.48 \\ 
##    \hline
## \end{tabular}
## \end{table}

table.usual<-data.frame(alpha.opt, rmse.usual, forecast.vol, VaR.dd[length(p),], ES.dd[length(p),])
table.elastic<-data.frame(omega.opt, rmse.elastic.opt, forecast.vol.elastic, VaR.dd.elastic[length(p), ], ES.dd.elastic[length(p),])

model.risk.VaR <-matrix(0, ncol=number, nrow=length(p))
model.risk.ES <-matrix(0, ncol=number, nrow=length(p))
for(num in 1:number){
  VaR<-data.frame(VaR_rateG[, num], VaR_ratetG1[, num], VaR.dd[, num], VaR.dd.res[, num])
  ES<-data.frame(ES_rateG[, num], ES_ratetG1[, num], ES.dd[, num], ES.dd.res[, num])
  model.risk.VaR[, num]<-apply(VaR, 1, max)/apply(VaR, 1, min)
  model.risk.ES[, num]<-apply(ES, 1, max)/apply(ES, 1, min)
}

model.risk.VaR.elastic<-matrix(0, ncol=number, nrow=length(p))
model.risk.ES.elastic <-matrix(0, ncol=number, nrow=length(p))
for(num in 1:number){
  VaR<-data.frame(VaR_rateG[, num], VaR_ratetG1[, num], VaR.dd.elastic[, num], VaR.dd.elastic.res[, num])
  ES<-data.frame(ES_rateG[, num], ES_ratetG1[, num], ES.dd.elastic[, num], ES.dd.elastic.res[, num])
  model.risk.VaR.elastic[, num]<-apply(VaR, 1, max)/apply(VaR, 1, min)
  model.risk.ES.elastic[, num]<-apply(ES, 1, max)/apply(ES, 1, min)
}

MR.usual<-data.frame(model.risk.VaR[length(p),], model.risk.ES[length(p),])
colnames(MR.usual)<-c("VaR", "ES")
rownames(MR.usual)<-main.names
MR.usual

##            VaR       ES
## SP500 2.186598 2.929684
## IBM   1.831004 2.292130
## MSFT  2.078809 2.670553
## AAPL  2.378313 2.964999
## GOOG  1.931144 2.507516

MR.elastic<-data.frame(model.risk.VaR.elastic[length(p),], model.risk.ES.elastic[ length(p),])
colnames(MR.elastic)<-c("VaR", "ES")
rownames(MR.elastic)<-main.names
MR.elastic

##            VaR       ES
## SP500 2.072867 2.792636
## IBM   1.753915 2.195627
## MSFT  1.967009 2.526928
## AAPL  2.296205 2.862637
## GOOG  1.844097 2.395500

MR<-data.frame(MR.usual, MR.elastic)
xtable(MR)

## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
##   \hline
##  & VaR & ES & VaR.1 & ES.1 \\ 
##   \hline
## SP500 & 2.19 & 2.93 & 2.07 & 2.79 \\ 
##   IBM & 1.83 & 2.29 & 1.75 & 2.20 \\ 
##   MSFT & 2.08 & 2.67 & 1.97 & 2.53 \\ 
##   AAPL & 2.38 & 2.96 & 2.30 & 2.86 \\ 
##   GOOG & 1.93 & 2.51 & 1.84 & 2.40 \\ 
##    \hline
## \end{tabular}
## \end{table}

temp1 <- data.frame (p, as.data.frame(model.risk.VaR))
colnames (temp1) <- c ("p","S&P 500","IBM","MSFT","AAPL", "GOOG")
temp1 <- gather (temp1, "Stock", "Model_risk", 2:6)
temp2 <- data.frame (p, as.data.frame(model.risk.ES))
colnames (temp2) <- c ("p","S&P 500","IBM","MSFT","AAPL", "GOOG")
temp2 <- gather (temp2, "Stock", "Model_risk", 2:6)
plot1 <- ggplot(temp1) +
  geom_point (aes (x = p, y = Model_risk, color = Stock, shape = Stock)) +
   ylab ("Model risk of VaR: no penalty")+ylim(c(1.5, 4))
plot2 <- ggplot(temp2) +
  geom_point (aes (x = p, y = Model_risk, color = Stock, shape = Stock)) +
   ylab ("Model risk of ES: no penalty")+ylim(c(1.5, 5))
grid.arrange(plot1, plot2, ncol=2)

temp1 <- data.frame (p,as.data.frame(model.risk.VaR.elastic))
colnames (temp1) <- c ("p","S&P 500","IBM","MSFT","AAPL", "GOOG")
temp1 <- gather (temp1, "Stock", "Model_risk", 2:6)
temp2 <- data.frame (p, as.data.frame(model.risk.ES.elastic))
colnames (temp2) <- c ("p","S&P 500","IBM","MSFT","AAPL", "GOOG")
temp2 <- gather (temp2, "Stock", "Model_risk", 2:6)
plot1 <- ggplot(temp1) +
  geom_point (aes (x = p, y = Model_risk, color = Stock, shape = Stock)) +
  ylab ("Model risk of VaR: elastic net")+ylim(c(1.5, 4))
plot2 <- ggplot(temp2) +
  geom_point (aes (x = p, y = Model_risk, color = Stock, shape = Stock)) +
  ylab ("Model risk of ES: elastic net")+ylim(c(1.5, 5))
grid.arrange(plot1, plot2, ncol=2)

model.risk.V<-list()
par(mfrow=c(1, 5))
for(j in 1:5){
  model.risk.V[[j]]<-data.frame(model.risk.VaR.elastic[,j], model.risk.VaR[, j])
  colnames(model.risk.V[[j]])<-c("REG", "NR")
  boxplot(model.risk.V[[j]], col=c(2, 3), main=main.names[j])
}

model.risk.E<-list()
par(mfrow=c(1, 5))
for(j in 1:5){
  model.risk.E[[j]]<-data.frame(model.risk.ES.elastic[,j ], model.risk.ES[, j])
  colnames(model.risk.E[[j]])<-c("REG", "NG")
  boxplot(model.risk.E[[j]], col=c("yellow", 5), main=main.names[j])
}

model.risk.0.01.nonPen<-data.frame(model.risk.VaR[length(p),], model.risk.ES[ length(p),])
model.risk.0.01<-data.frame(model.risk.VaR.elastic[length(p),], model.risk.ES.elastic[ length(p),], model.risk.0.01.nonPen)
rownames(model.risk.0.01)<-main.names
colnames(model.risk.0.01)<-c("VaR(REG)", "ES(REG)", "VaR", "ES")

#library(vioplot)
par(mfrow=c(1, 5))
for(i in 1:number){
  vioplot(model.risk.V[[i]], col=c(2, 3), names=c("REG", "NR"))
  title(paste(main.names[i], "(VaR)"), ylab="model risk")
}

par(mfrow=c(1, 5))
for(i in 1:number){
  vioplot(model.risk.E[[i]], col=c("yellow", 5), names=c("REG", "NR"))
  title(paste(main.names[i], "(ES)"), ylab="model risk")
}

EVT for VaR and ES forecasting: Reg:

res<-tail(res, 1000)
res.elastic<-tail(res.elastic, 1000)
m<-seq(150, 300, 10)
m.opt<-rep(0, number)
for(j in 1:number){
  res.use<-res[, num]
  mse<-rep(0, length(m))
  for(i in 1:length(m)){
    res.order<-sort(res.use)[1:m[i]]
    y<-log((-1)*res.order)
    x<-log(seq(1, m[i], 1)/1000)
    fit<-lm(y~x)
    mse[i]<-mean(fit$residuals^2)
  }
  m.opt[j]<-m[which.min(mse)]
}
m.opt

## [1] 150 150 150 150 150

m.opt.elastic<-rep(0, number)
for(j in 1:number){
  res.use<-res.elastic[, num]
  mse<-rep(0, length(m))
  for(i in 1:length(m)){
    res.order<-sort(res.use)[1:m[i]]
    y<-log((-1)*res.order)
    x<-log(seq(1, m[i], 1)/1000)
    fit<-lm(y~x)
    mse[i]<-mean(fit$residuals^2)
  }
  m.opt.elastic[j]<-m[which.min(mse)]
}
m.opt.elastic

## [1] 150 150 150 150 150

a.reg.inv<-rep(0, number)
a.reg.inv.elastic<-rep(0, number)
for(j in 1:number){
  res.use<-sort(res[, j])[1:m.opt[j]]
  y<-log((-1)*res.use)
  x<-log(seq(1, m.opt[j], 1)/1000)
  fit<-lm(y~x)
  a.reg.inv[j]<-(-1)*fit$coefficients[2]
}

for(j in 1:number){
  res.use<-sort(res.elastic[, j])[1:m.opt.elastic[j]]
  y<-log((-1)*res.use)
  x<-log(seq(1, m.opt.elastic[j], 1)/1000)
  fit<-lm(y~x)
  a.reg.inv.elastic[j]<-(-1)*fit$coefficients[2]
}
a.reg.inv

## [1] 0.4856916 0.4728223 0.4574394 0.3882623 0.4347708

a.reg.inv.elastic

## [1] 0.4673441 0.4725276 0.4498897 0.3823091 0.4295732

p<-0.1
VaR.reg<-rep(0, number)
VaR.reg.elastic<-rep(0, number)
for(j in 1:number){
  VaR.reg[j]<-(-1000)*sort(res[, j])[100]*forecast.vol[j]*(0.1/0.01)^a.reg.inv[j]
  VaR.reg.elastic[j]<-(-1000)*sort(res.elastic[, j])[100]*
                   forecast.vol.elastic[j]*(0.1/0.01)^a.reg.inv.elastic[j]
}
VaR.reg; VaR.reg.elastic

## [1] 47.24521 65.77080 75.12274 83.91127 69.36910

## [1] 42.80098 59.26378 65.42395 76.52281 62.27243

ES.reg<-VaR.reg*(1/(1-a.reg.inv))
ES.reg.elastic<-VaR.reg.elastic*(1/(1-a.reg.inv.elastic))

Hill:

nc<-seq(5, 300, 5)
hill.est<-function(y, nc){
  y.use<-sort(y)[1:nc]
  a.hill<-nc/sum(log(y.use/y.use[nc]))
  return(a.hill)
}
for(j in 1:number){
  a.hill<-rep(0, length(nc))
  for(i in 1:length(nc)){
    a.hill[i]<-hill.est(res[, j], nc[i])
  }
  plot(nc, a.hill, type="b", main=main.names[j])
}

nc.opt<-100
for(j in 1:number){
  a.hill<-rep(0, length(nc))
  for(i in 1:length(nc)){
    a.hill[i]<-hill.est(res.elastic[, j], nc[i])
  }
  plot(nc, a.hill, type="b", main=main.names[j])
}

nc.opt.elastic<-100
a.hill.inv<-rep(0, number)
a.hill.inv.elastic<-rep(0, number)
for(j in 1:number){
  a.hill.inv[j]<-1/hill.est(res[, j], nc.opt)
  a.hill.inv.elastic[j]<-1/hill.est(res[, j], nc.opt.elastic)
}
VaR.hill<-rep(0, number)
VaR.hill.elastic<-rep(0, number)
for(j in 1:number){
  VaR.hill[j]<-(-1000)*sort(res[, j])[100]*forecast.vol[j]*(0.1/0.01)^a.hill.inv[j]
  VaR.hill.elastic[j]<-(-1000)*sort(res.elastic[, j])[100]*
                   forecast.vol.elastic[j]*(0.1/0.01)^a.hill.inv.elastic[j]
}
ES.hill<-VaR.hill*(1/(1-a.hill.inv))
ES.hill.elastic<-VaR.hill.elastic*(1/(1-a.hill.inv.elastic))
ES.hill; ES.hill.elastic

## [1]  94.32527 101.66041 126.63569 166.99158 136.43301

## [1]  89.13979  91.66483 112.22022 154.38967 123.95010

EVT.NP<-data.frame(VaR.hill, VaR.reg, ES.hill, ES.reg)
EVT.REG<-data.frame(VaR.hill.elastic, VaR.reg.elastic, ES.hill.elastic, ES.reg.elastic)
EVT.NP; EVT.REG

##   VaR.hill  VaR.reg   ES.hill    ES.reg
## 1 47.92611 47.24521  94.32527  91.86162
## 2 58.65128 65.77080 101.66041 124.76021
## 3 71.45808 75.12274 126.63569 138.45962
## 4 93.95547 83.91127 166.99158 137.16870
## 5 73.60442 69.36910 136.43301 122.72739

##   VaR.hill.elastic VaR.reg.elastic ES.hill.elastic ES.reg.elastic
## 1         45.29139        42.80098        89.13979        80.3539
## 2         52.88450        59.26378        91.66483       112.3543
## 3         63.32371        65.42395       112.22022       118.9288
## 4         86.86519        76.52281       154.38967       123.8853
## 5         66.87000        62.27243       123.95010       109.1681

EVT<-data.frame(EVT.NP, EVT.REG)
rownames(EVT)<-main.names
EVT

##       VaR.hill  VaR.reg   ES.hill    ES.reg VaR.hill.elastic VaR.reg.elastic
## SP500 47.92611 47.24521  94.32527  91.86162         45.29139        42.80098
## IBM   58.65128 65.77080 101.66041 124.76021         52.88450        59.26378
## MSFT  71.45808 75.12274 126.63569 138.45962         63.32371        65.42395
## AAPL  93.95547 83.91127 166.99158 137.16870         86.86519        76.52281
## GOOG  73.60442 69.36910 136.43301 122.72739         66.87000        62.27243
##       ES.hill.elastic ES.reg.elastic
## SP500        89.13979        80.3539
## IBM          91.66483       112.3543
## MSFT        112.22022       118.9288
## AAPL        154.38967       123.8853
## GOOG        123.95010       109.1681